function vmd_images(fbase, fnum, OmegaBounce)
%
%   vmd_images(fbase, fnum, OmegaBounce)
%
%   Reads file with density distribution and file with x-velocity.
%   Creates bunch of images.
%
%   Input:
%        fbase   base of file number (runname)
%        fnum    file number
%
%   Domain information and colormap limits are edited in the file.
%   (OmegaBounce = 5.3 at t=200 in the test run)


%---------------------------------------------------------

% system size

jx=48; jy=1; jz=16; jv=256; jdir=2;

% physics paramters

vmax=10; Lx=450; kz=0.34; nz=1; vPhase=3.53235; u0=1.24; phi00=0.01;
SpeckleWidth = 99.91139;

% colormap bounds

rhomax = 0.02;  uxmax = 0.005;  uzmax = 0.075;
phimax = 0.2;   exmax = 0.006;  ezmax = 0.06;
fzvmax = 0.02;

%---------------------------------------------------------

dv = 2*vmax/jv;
dx = Lx/(jx-1);       % non-periodic BC
dz = 2*pi*nz/(jz*kz);

dkz = kz/nz;
dkx = 2*pi / Lx;

v  = (0:jv-1)*dv - vmax - vPhase;
z  = (0:jz-1)*dz;
x  = (0:jx-1)*dx - Lx/2;

f0 = ( 1/sqrt(2*pi) )*exp( - 0.5*(v + vPhase).**2 );

xmod = exp( -20*(x/SpeckleWidth).**2 );
          
phi0=zeros(jz,jx); 
for iz=1:jz
  phi0(iz,:) = xmod * phi00 *exp(i*kz*z(iz)) ;
end

phi0FT = fft2(phi0);


%phi = real(phi0FT);
%disp([ min(min(phi)), max(max(phi))])
%imagesc(flipud(phi))
%return

%---------------------------------------------------------

%-- reading density distribution

  fn = num2str(fnum,'%04d');

  fname = [fbase, '.dfu.', fn];

  fid = fopen(fname, 'rb');
  dfu = fread(fid, jv*jx*jz*jdir, 'double');
  fclose(fid);

  dfu = reshape(dfu, jz,jv,jx,jdir);


%-- distrubution over velocities --

  ix = jx/2+1;
  fzv =  dfu(:,:,ix, 1) + dfu(:,:,ix, 2);

  fzv =  reshape(fzv, jz, jv)';
  
  fzv =  fzv + meshgrid(f0, 1:jz)';

  vlimit = sqrt(2*OmegaBounce /SpeckleWidth/0.37233 /kz);

 
  indv = find( abs(v) < vlimit );

  fzv = fzv(indv, :);

  printf("trapped limits :        [%8.4f, %8.4f]\n",  min(min(fzv(:,:))), max(max(fzv(:,:))) )


  figure(7, "Name", "Trapped")
  imagesc(flipud(fzv))
  one_image(flipud(fzv), 0, fzvmax, [fbase, '.fvz.', fn, '.png' ])



%-- integrate over "v" to find density

  rho = sum(sum(dfu, 2), 4); 
  rho = reshape(rho, jz, jx)*dv;

  printf("density limits :        [%8.4f, %8.4f]\n",  min(min(rho(:,:))), max(max(rho(:,:))) )

  figure (1, "Name", "Density")
  imagesc(flipud(rho));
  one_image(flipud(rho), -rhomax, rhomax, [fbase, '.rho.', fn, '.png' ])


%-- integrate over "v" to find z-velocity

  for iv=1:jv
    dfu(:,iv,:,:) = dfu(:,iv,:,:) + 0.5*f0(iv);
    uz (:,iv,:,:) = dfu(:,iv,:,:).*v(iv);
  end

  uz = sum(uz,2)./sum(dfu,2);
  uz = reshape(uz, jz, jx,jdir);

  uz1 = uz(:,:,1) + vPhase;
  uz2 = uz(:,:,1) + vPhase;

  clear dfu

  printf("z-velocity (1) limits : [%8.4f, %8.4f]\n", min(min(uz1(:,:))), max(max(uz1(:,:))) )
  printf("z-velocity (2) limits : [%8.4f, %8.4f]\n", min(min(uz2(:,:))), max(max(uz2(:,:))) )

  figure (2, "Name", "z-Velocity (fluid 1)")
  imagesc(flipud(uz1));
  one_image(flipud(uz1), -uzmax, uzmax, [fbase, '.Uz.', fn, '.png' ])

  clear uz
  clear uz1
  clear uz2

%-- reading x-velocity

  fname = [fbase, '.ux.', fn];

  fid = fopen(fname, 'rb');
  ux = fread(fid, jx*jz*jdir, 'double');
  fclose(fid);

  ux = reshape(ux,jx,jz,jdir);

  ux1 = (ux(:,:,1) - u0)';
  ux2 = (ux(:,:,2) + u0)';

  printf("x-velocity (1) limits : [%8.4f, %8.4f]\n",  min(min(ux1(:,:))), max(max(ux1(:,:))) )
  printf("x-velocity (2) limits : [%8.4f, %8.4f]\n",  min(min(ux2(:,:))), max(max(ux2(:,:))) )

  figure(3, "Name", "x-Velocity (fluid 1)")
  imagesc(flipud(ux1));
  one_image(flipud(ux1), -uxmax, uxmax, [fbase, '.Ux.', fn, '.png' ])

  clear ux
  clear ux1
  clear ux2

%-- solve poisson equation 

  rho = fft2(rho);

  kx = [0:jx/2,-jx/2+1:-1] *dkx;
  kz = [0:jz/2,-jz/2+1:-1] *dkz;

  [kkx, kkz] = meshgrid(kx, kz);

  kk = kkx.*kkx + kkz.*kkz;

  phi = rho./kk;


  phi(:,1) = 0;
  phi = phi + phi0FT;
  phi(1,:) = 0;


  Ex  = -i*kkx.*phi;  
  Ez  = -i*kkz.*phi;  

  Ex  = real(ifft2(Ex));
  Ez  = real(ifft2(Ez));
  phi = real(ifft2(phi));

  printf("E_x limits :            [%8.4f, %8.4f]\n",  min(min(Ex(:,:))), max(max(Ex(:,:))) )
  printf("E_z limits :            [%8.4f, %8.4f]\n",  min(min(Ez(:,:))), max(max(Ez(:,:))) )
  printf("phi limits :            [%8.4f, %8.4f]\n",  min(min(phi(:,:))), max(max(phi(:,:))) )


  figure(4, "Name", "x-Efield")
  imagesc(flipud(Ex))
  one_image(flipud(Ex), -exmax, exmax, [fbase, '.Ex.', fn, '.png' ])

  figure(5, "Name", "z-Efield")
  imagesc(flipud(Ez))
  one_image(flipud(Ez), -ezmax, ezmax, [fbase, '.Ez.', fn, '.png' ])

  figure(6, "Name", "Potential")
  imagesc(flipud(phi))
  one_image(flipud(phi), -phimax, phimax, [fbase, '.phi.', fn, '.png' ])


end


%---------------------------------------------------------

function one_image(a, amin, amax, fname)

  a = max(a,amin);
  a = min(a,amax); 

  map = colormap(jet(256));

  sc   =  uint8(255*(a-amin)/(amax-amin)) + 1;
  imwrite(sc, map, fname);

end


%---------------------------------------------------------
